Ultralow Interfacial Thermal Resistance of Graphene Thermal Interface Materials with Surface Metal Liquefaction

Highlights A three-tiered thermal interface materials was proposed with the through-plane thermal conductivity up to176 W m−1 K−1 and contact thermal resistance as low as 4–6 K mm2 W−1 (double sides). The liquid metal acts as a buffer layer to connect vertically aligned graphene with the rough heater/heat sink, improving effective contact thermal conductance by more than an order of magnitude. Supplementary Information The online version contains supplementary material available at 10.1007/s40820-022-00979-2.


Supplementary
Comparison of intrinsic through-plane conductivity between LM-VAGM and reported thermally conductive materials based on vertically aligned structure

Flash Method
The thermal conductivity of LM-VAGM was measured by a laser flash method using LFA 467 HyperFlash ® system (NETZSCH, Germany). To carry out the measurement, the samples were cut into a square plate with the size of 10 × 10 × 0.8 mm 3 . All samples are coated by a graphite layer (≈ 5 μm) on both upper and lower sides.
In the actual operation, one side of the sample set in the holder was transiently heated by a xenon lamp pulse with a pulse width of 10 μs, and an infrared radiation detector recorded the rise of surface temperature at another side. The thermal diffusivity was evaluated by analyzing the fitted curve of temperature evolution versus sampling time based on Eq. S1: where α, d, and t50 are the calculated thermal diffusivity, sample thickness, and the half diffusion time, respectively. The α was obtained by measuring three samples separately three times. Accordingly, the thermal conductivity (κ) could be calculated using Eq. S2: where Cp is the specific heat capacity and is measured using LFA 467 with pure graphite as reference (25 °C); ρ is the density of the samples and is calculated from the ratio of mass to volume. The Cp and the ρ were obtained by taking the average value of three separate measurements.

S2.1 Uncertainty Analysis of LM-VAGM Thermal Conductivity
Based on the measurement principle of the flash method and Eqs. S1 and S2, the uncertainty of the LM-VAGM thermal conductivity mainly comes from the following four factors: μ1uncertainty calculated from repeated measurements of sample thermal diffusivity (αLM-VAGM) by statistical analysis (A type). μ2uncertainty introduced by measuring the thickness of a sample (B type). μ3uncertainty introduced by measuring the specific heat capacity of a sample (B type). μ4uncertainty introduced by measuring the density of a sample (B type). Based on the result in Table S2, the standard deviation of αLM-VAGM can be calculated using Eq. S3:

S2.2 Calculation of μ1
Accordingly, the uncertainty calculated from repeated measurements of sample thermal diffusivity can be obtained using Eq. S4:

S2.3 Calculation of μ2
The sample thickness was measured by a micrometer and the calibration certificate indicates an uncertainty of 1.3 μm (k = 2). Based on the sample thickness of 800 μm, the relative standard uncertainty (μthickness) can be obtained: Based on the Eq. S1, uncertainty introduced by measuring the thickness of a sample can be obtained by the following:

S2.4 Calculation of μ3
Based on the instruction of the instrument, the test accuracy of the specific heat capacity is ±5%.
Considering the uncertainty according to the normal distribution, the confidence probability is 99.73%. Therefore, the uncertainty introduced by measuring the specific heat capacity can be obtained:

S2.5 Calculation of μ4
The density of the sample is calculated from the ratio of mass (m) to volume (V). The sample mass was measured using analytical balance, with an uncertainty of 0.1 mg (k = 2). Based on the sample mass of 53.6 mg, the relative standard uncertainty (μm) can be obtained: The width and length were measured by a micrometer, with an uncertainty of 1.3 μm (k = 2). Based on the sample width/length of 10 mm, the relative standard uncertainty (μwidth/length) can be obtained: Accordingly, the uncertainty introduced by measuring the density can be obtained:

S2.7 Calculation of Expanded Uncertainty: Urel
Taking the coverage factor k = 2, Urel can be obtained: The Rtotal for a TIM was measured using the modified ASTM D5470 method, as schematically illustrated in Fig. S2. BLT (bond line thickness) is the thickness of test TIM under compression. X1 -X6 is the location of the temperature measurement point on the heat flux meter bar. For the Longwin 9389 equipment, X1 = X2 = X 4 = X 5 = X 6 = 25 mm; X 3 = X 4 = 15 mm. Tc1 -Tc3 and Th1 -Th3 are the measured temperature using six resistance temperature detectors (RIDs). Tc and Th are the calculated temperature of the upper and lower surfaces of the test TIM, respectively. Qc and Qh are the heat flux of the heat meter bar and cold meter bar, respectively. Qc and Qh can be extracted from the one-dimensional heat equation, as shown in Eqs. S13 and S14: where κCu and ACu are the thermal conductivity and the heat transfer area of the heat flux meter bar, which is made of pure copper, plated with a nanometer-thick layer of nickel. The average heat flux can be calculated using Eq. S15: When the test reaches a steady state, Eqs. S16 and S17 can be used to describe the relationship between the temperature and the location of the temperature measurement point.
Accordingly, the upper (Tc) and lower surfaces (Th) of the test TIM can be obtained based on Eqs. S18 and S19, and the Rtotal can be calculated using Eq. S20, where ATIM is the area of the tested TIMs.
Based on the above principle, we selected three sets of samples with the different original thicknesses (200, 500, and 800 μm) and measured the Rtotal from 10 to 100 psi, as the results shown in Fig. S3a, b, respectively. To obtain the ⊥ under compression and the corresponding Rcontact, we perform a linear fit to the measured Rtotal versus the BLT at each pressure. Fig. S3c, d exhibited the typical results for VAGM and LM-VAGM at 30 psi, in which the inverse of the calculated slope is ⊥ and the intercept is Rcontact. Based on the same data processing method, we got the ⊥ and Rcontact of VAGM and LM-VAGM from 10 to 100 psi.

S3.1 Uncertainty Analysis of the total Resistance of the LM-VAGM
Based on the measurement principle of the ASTM D5470 method using the test equipment of LW 9389 (Long Win), the uncertainty of the total thermal resistance for the LM-VAGM mainly comes from the following three factors: μ5uncertainty calculated from repeated measurements of sample total thermal resistance (RLM-VAGM) by statistical analysis (A type).

μ6uncertainty introduced by the random error of test equipment (B type).
μ7uncertainty introduced by measuring the lateral size of a sample (B type). Based on the result in Table S3, the standard deviation of RLM-VAGM can be calculated using Eq. S21:

S3.2 Calculation of μ5
Accordingly, the uncertainty calculated from repeated measurements of sample thermal conductivity can be obtained:

S3.3 Calculation of μ6
The random error of the test equipment mainly comes from three aspects, including temperature difference measurement error, contact condition error, and instrument accumulated error. According to the device manual (Long Win LW-9389), the values of the three random errors are ±12.0% (temperature difference measurement error), ±3.0% (contact condition error), and ±1.0% (instrument accumulated error), respectively. Treating all these errors as mean distributions, the uncertainty introduced by the random error of test equipment can be obtained:

S3.4 Calculation of μ7
The lateral size (width × length) was measured by a micrometer and the calibration certificate indicates an uncertainty of 1.3 μm (k = 2). Based on the sample width/length of 25.4 mm, the relative standard uncertainty (μthickness) can be obtained: Accordingly, the uncertainty introduced by measuring the thickness of a sample can be obtained by the following: 7 = √( width/length ) 2 + ( width/length ) 2 = 0.037% (S25)

Multi-Block Lattice Boltzmann Method
In order to in-depth elaborate on the interfacial heat transfer process, we used the multi-block lattice Boltzmann method to calculate the temperature distribution of the mating interface formed by a TIM and a heat sink (Fig. S4a). The detailed model setup and solution steps are as follows.

S4.1 Rough Surface Generation
The numerical model is based on the practical topography of the specimens. Firstly, the high-resolution microscope was used to measure the surface topography of the TIM and the heat sink, as the results shown in Fig. S4b, c, respectively. The measured area is 5 mm × 5 mm and the unit of height is μm. Then, the measured surface topography data was used to generate the lower surface of the TIM and the upper surface of the heat sink.
It should be noted that the measured data is quite massive that could not be inputted into the software directly, so 1/64 data volume was used to reconstruct rough surfaces. Here, 1/64 means that the used data volume accounts for 1/64 of the original measured data volume. Accordingly, the numerical model of the TIM/heat sink mating interface could be constructed by importing the surface topography data into a pre-processing software named ANSA via the "points→coons surfaces→volume→mesh" flow path. As a result, the coons surface generation was realized through the self-developed Python code read in the ANSA script interface. Figure S4d, e exhibits the imported data points and generated coons surfaces. The distance of every point in both x and y directions is 0.031807 mm. Here, we used 24964 points to generate a rough surface, and the area of every generated rough surface was 4.998 × 4.998 mm 2 . It can be seen that every small coons surface was generated by four adjacent points, and all small coons surfaces formed an entire continuous rough surface through the "Topo function" of the ANSA software. Fig. S4 (a) Schematic illustrating the VAGM as a TIM in contact with the rough surface of the heat S11/S14 sink. The measured surface roughness of (c) VAGM and (d) the heat sink. (d) The imported data points of the VAGM based on the measured surface roughness, with the reconstructed surface of the TIM shown in (e).

S4.2 Mesh generation
Based on the reconstructed rough surfaces of the TIM and the heat sink, the numerical model of the TIM/heat sink mating interface could be built, as shown in Fig. S5a. The height of the TIM and the heat sink are 0.8 mm and 2 mm, and the lateral size of them is 4.998 × 4.998 mm 2 . When generating meshes, the mating interfaces were meshed firstly, and the imported points were treated as nodes of a single quadrilateral mesh. Then, the quadrilateral meshes were mapped into structured hexahedral meshes. Finally, the whole model has 1469520 structured elements, as the corresponding mating interface shown in Fig. S5b.

S4.3 Materials property settings
The simulation is finished by finite element software ABAQUS. According to the experimental data, the material of the heat sink is set to copper (Cu) with the thermal conductivity of 398 W m -1 K -1 and the compression modulus of 110 GPa. For the case without liquid metal between TIM and heat sink, we set TIM material properties based on experimental data from VAGM (thermal conductivity: 210 W m -1 K -1 ; compression modulus: 2.25 MPa). For the case with liquid metal as a buffer layer, we did not directly use the experimental data of LM-VAGM, but reduced the compressive modulus by 10% based on the material properties of VAGM.

S4.4 Boundary conditions and solution
After setting the material properties, the contact computation between the lower surface of the TIM and the upper surface of the heat sink was conducted. In this stage, the thermal expansion between S12/S14 them was neglected, and different pressures (10, 20, and 30 psi) were loaded on the top surface of TIM. When the contact computation was finished, the actual contact area of the mating interface under different pressure can be obtained. Then, a steady-state heat analysis was conducted on the deformed model. In this stage, a heat flux of 12.4 W cm -2 was loaded at the top surface of TIM, and the bottom surface was kept at 68 ℃. When the computation converges, we can obtain the temperature distribution of the mating interface, as the results shown in Fig. S5c, d, respectively. In order to in-depth study the TIM performance in an integrated system using LM-VAGM, a commercial computational fluid dynamics software (Icepak) was employed to simulate the heat transfer process in an electronic system with the power of the device (heater) of 200 W. The model implementation is presented in Fig. S6a, and the detailed settings of the heater and the heat sink can be found in Table S4. When the simulation started, the background temperature was set to 20 °C at 1 atm in the atmosphere, and the entrance temperature of cooling water was 20 °C with a volume flow of 1.26 cm 3 s -1 . Fig. S6b-d shows the steady-state temperature profiles of the simulated system with TIMs, demonstrating a better heat dissipation capability of our LM-VAGM compared to that with VAGM or Carbonaut thermal pad as the TIM. In addition, based on the equation: Rcontact = Rtotal -Rbulk = BLT/κeff -BLT/κbulk, the total (Rtotal) and contact thermal resistances (Rcontact) of the three TIMs were calculated, with the detailed parameters and results shown in Table S5.